Stability of K-montmorillonite hydrates: Hybrid MC simulations 
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NPzzT and jiPzzT simulations of K-montmorillonite hydrates were performed employing hybrid 
Monte Carlo simulations. Two condition sets were studied, P=l atm and T= 300 K (ground level 
conditions), and P=600 atm and T= 394 K; this last condition mimics a burial depth close to 4 
km. For these conditions, swelling curves as a function of the reservoir water vapor pressure were 
built. We found the single layer K-montmorillonite hydrate stable for high vapor pressures for both, 
burial and ground level conditions. A simple explanation for this high stability is given. 

PACS numbers: 



I. INTRODUCTION 

Clays are layer type aluminosilicate minerals, existent 
everywhere in nature and industry, hence the importance 
of a detailed understanding of their physics and chem- 
istry. They are used as building materials, ceramics, and 
catalysts; they are employed in cosmetics, as rheological 
modifiers for paints, and in technological processes such 
us oil well drilling. In this last application, the control of 
their stability is a key for drilling success. During it, the 
use of water based muds induce destabilization of shale 
and clay formations that would disintegrate, or heave, 
upon contact with water. 

One way to maintain stability of shales during the 
drilling process is by addition of potassium salts to 
drilling muds. This helps to avoid fluid loss and wa- 
ter infiltration. These kind of muds, that contain potas- 
sium ions dissolved in the water phase, are widely used 
for drilling water-sensitive shales, specially hard, brittle 
shales. Potassium cations in these systems replace ions 
such as sodium found in most shales to produce less hy- 
drated clays with significantly reduced swelling potential. 
These ions also help to hold the cuttings together, min- 
imizing their dispersion into finer particles. From all of 
the previous facts, a good understanding of the role of 
potassium in the swelling clays such as montmorillonite 
is mandatory, in particular at basin conditions of hard 
experimental implementation. 

Besides the experimental studies, computer simula- 
tions are essential components of research on clay- water- 
cation systema-^i^i^i'^i^i^i'''i^i^i-^°i-^-^i-^^i^^. These simulations 
give microscopic insights that are difficult to access ex- 
perimentally. There are, however, not many computa- 
tional studies on the swelling of montmorillonite hydrates 
for potassium interlayer cation s rfi^^-^^'^''''^^'^^'^" . In ad- 
dition, excluding Hensen et ali^ and Tambach et alm^ 
works, these papers report solely NP^zT and NVT sim- 
ulations, where the number of water molecules per in- 
terlaminar space is someway arbitrarily fixed, and hence, 
they do not show the whole picture of swelling. A good 
example of this is that hysteresis is naturally predicted by 
sampling in an open ensemble^i2i2i, without the need of 
measuring properties such as water chemical potential^i, 
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immersion energiesiii, or swelling free energie3«ia*. Fi- 
nally, this paper focuses on the stability of the different 
hydrates in contact with several reservoirs, which differ 
in temperature, pressure, and water activity. 

The importance of potassium as swelling inhibitor of 
clays and the above mentioned reasons motivated us to 
study the microscopic mechanisms underlying the behav- 
ior of K-montmorillonite hydrates at equilibrium with 
different reservoirs. We performed simulations of these 
systems in the NP^zT and jJiPzzT ensembles, considering 
explicitly two clay layers in the simulation box to avoid 
finite size effects. The simulations were carried out for 
two condition sets. One at ground level, with P= 1 atm 
and T= 298K, and the other one with P= 600 atm and 
T= 394K, which corresponds to an average burial depth 
close to 4 km. 

The paper is organized as follows. In Sec.^ we briefly 
describe the models and the methodology employed for 
performing the simulations. The results are shown in 
Sec. mil Finally, Sec. IIVI discusses the main results and 
extracts some conclusions. 



II. METHODOLOGY 

A. The model 

A 4 X 2 layer of Wyoming type montmorillonite clay 
was built up by replication of the unit cell given by Skip- 
per et alM. This layer has = 21.12 A, Ly ^ 18.28 
A and Lz — 6.56 A dimensions. The Wyoming type 
montmorillonite was obtained by isomorphous substitu- 
tions of trivalent Al atoms of the octahedral sites by 
divalent Mg atoms, and tetravalent Si by trivalent Al 
atoms. The unit cell formula of this clay is given by 
Ko.75nH20(Si7.75Alo.25)(Al3.5Mgo.5)02o(OH)A Size ef- 
fects were avoided by considering two layers in the simu- 
lation box**. Periodic boundary conditions were imposed 
on the three space directions. The initial configuration 
consists of water molecules randomly placed in the in- 
terlaminar spaces, and six potassium ions distributed in 
the interlayer midplanes. These counterions balance the 
negative charge of the clay framework keeping the system 
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TABLE I: Lennard- Jones parameters for H20-clay-K inter- 
actions. Exceptions are the K-O and K-H interactions. 



Sites 



e (kcal/mol) 



^ (A) 



O 
H 
K 

Si 
Al 
Ms 



0.1550 
0.0000 
3.630 
3.150 
3.150 
3.150 



3.1536 
0.0000 
2.4500 
1.840 
1.840 
1.840 



electroneutral. 

The rigid TIP4P model was used for water 
moleculesSSi^ and the water clay interactions were taken 
from Boek et al^ Here, site to site intermolecular in- 
teractions are given by electrostatic and Lennard-Jones 
contributions, 
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where subindexes i and j are for molecules, and a and 
h run over all sites of each molecule, qa and qi, are the 
corresponding site charges, eab and Gab are site to site 
specific Lennard-Jones parameters and rab is the inter- 
site distance. The Lennard-Jones parameters for single 
sites are shown in Table Here, those parameters for Si 
were taken from Marry et alm^, and parameters for Al 
and Mg were assumed to be equal to those of Si. The 
site to site Lennard-Jones parameters are given by the 
Lorentz-Berthelot rules 
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On the other hand, the K-H2O interactions and those 
between the oxygens of the clay and potassium ions are 
based on the ones proposed by Bounds^^. Bounds poten- 
tial is chosen since, while simple, it produces K-TIP4P 
radial distribution functions in agreement with available 
experimental data and close to those obtained by hybrid 
quantum mechanics/ molecular mechanics (QM/MM) 
simulations, which naturally account for the many body 
contributions to the potentiaiSl. That is, the K-0 radial 
distribution function peaks at 2.86 A leading to a first 
shell oxygen coordination number of 7.6, while high ac- 
curacy QM/MM simulations performed at density func- 
tional theory level (LANL2DZ basis set) give 2.81 A 
of K-0 distance and 8.3 of coordination number. Ex- 
perimental results give K-0 distances between 2.7 and 
3.1 A and coordination numbers in the wide range of 
^_g28j29_ jjg fitted the following pair potential for the 
K-H2O dispersion- repulsion contribution obtained from 



ab-initio calculations, 

UK-H2O = Ako exp {-bKorxo) - CKo/r%o 

-DKo/rKo + ^KH exp {-bxHrKHx ) 
+AKHe^v{-bKHrKH2), (4) 

yielding Ako = 53884.0 kcal/mol, Bko = 

3.3390 A , 

Cko = 438.0 kcalA'^/mol, Dko =-638.0 kcalA^mol, 

Akh = 5747.0 kcal/mol and bxH = 3.4128 A~\ This 
intersite potential, although somewhat more complicated 
than the Lennard-Jones type, produces a much better 
match to the ab-initio data2&. This is clearly seen if one 
tries to fit equation 10} with a Lennard-Jones type po- 
tential. In fact, parameters shown in tabled for potas- 
sium were obtained by this way, yielding a not very good 
match although reproducing the depth and position of 
the pair-potential minimum. We observed that the dis- 
crepancies are very pronounced at short distances, where 
the Lennard-Jones potential shows a much harder be- 
havior. This explains why Boek et al. found a very 
large dehydrated interlaminar space when they employed 
a Lennard-Jones type pair potential for K-0^^. Natu- 
rally, they overcome this difficulty by employing the pair 
potential proposed by Bounds^^. 

Nevertheless, since it is crucial for the hybrid Monte 
Carlo simulations to keep the energy fluctuations as low 
as possible in order to enlarge the acceptation rata^, it is 
convenient to avoid employing relatively long range pair 
potential contributions such as ^ r~'^, if no Ewald treat- 
ment is applied on them. Hence, we refitted to equation 
0]the following expression 

UK-H20 = Axoexp {-bKoTRo) - CKo/r'ko 
+Akh exp {-bKHrKHi ) 
+Akh exp{-bKHrKH2), (5) 

by employing a Levenberg-Marquardt algorithm and con- 
sidering several K-H2O configurations. The procedure 
yields Ako = 120750.2 kcal/mol, bxo = 3.4110 A"\ 
Cko = 5153.8 kcalA^/mol, Akh = 2109.4 kcal/mol and 
bKH = 2.8515 A ^. We observed that avoiding the 
term the acceptance rate enlarges more than three times 
for a small (inner) time step of 0.8 fs. In general, both 
functions yield similar values of the K-H2O potential en- 
ergy. Minima are located practically at the same dis- 
tance although the depth of the fitted function is 5% 
larger. For larger distances this difference decreases. To 
check the obtained interaction potential, a NPT simula- 
tion containing 216 water molecules, a potassium cation, 
and a chloride anion was performed at P =1 atm and 
T=293 K. The corresponding K-0 and K-H radial distri- 
bution functions, g(r), and coordination numbers, n(r), 
were studied. These functions were observed to be very 
similar to those reported by Bounds^. That is, the K-0 
and K-H main peaks are located at 2.83 and 3.29 A re- 
spectively, which compare well with their corresponding 
values of 2.86 and 3.32 A^e. The coordination number for 
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the g(r) minimum was found at 7.7 A, in agreement with 
his value of 7.6 A^^. Moreover, all these values are even 
closer to the results obtained by QM/MM simulations^i. 
Hence, expression|Slseems to be suitable for our purposes. 

Finally, we should mention that electrostatic contribu- 
tions, ~ , were treated trough the implementation of 
the Ewald summation formalism. Here the convergence 
factor was fixed to 5.6/Lmin, where Lmm is the mini- 
mum simulation box side. There were set five reciprocal 
lattice vectors for the directions along the shortest sides 
and six vectors for the direction along the largest side^. 
The dispersion-repulsion contributions were corrected us- 
ing the standard methods for homogeneous fluidg.22, and 
a spherical cutoff of Lmm/2 was imposed. 

B. Simulations 

Simulations were performed employing the hybrid 
Monte Carlo (HMC) methodSLSa. This technique allows 
making global moves while keeping a high average ac- 
ceptance probability. Global moves are done as follows 
from molecular dynamics (MD), i. e., by assigning ve- 
locities and by using a particular scheme for integrat- 
ing the Newton's equations of motion. Velocities are as- 
signed randomly from a Gaussian distribution in corre- 
spondence with the imposed temperature, and in such 
a way that total momentum equals zero for both inter- 
laminar spaces. To fulfill detail balance condition, the 
discretization scheme must be time reversible and area 
preservingSa. In particular, we employed the multiple 
time scale algorithm given by Tuckerman and Bernc'^^. 
This algorithm has the property of splitting the forces 
into short and long range. The Lennard- Jones contribu- 
tion plus the real part of the electrostatic forces are set 
as short range and the reciprocal space contribution of 
the electrostatic forces is set as long range. To decrease 
time correlations a new configuration is generated each 
10 integration steps. The probability to accept this new 
configuration is given by 

P = min{l, exp(-/3AH)} (6) 

where ATi is the difference between the new and previ- 
ous configuration Hamiltonians, and (3 is the inverse of 
the thermal energy. The long time step is set to 8 times 
the short time step, and the short time step is chosen to 
obtain an average acceptance probability of 0.7'^''. This 
way, we obtained short time steps close to 1.0 and 0.5 
fs for systems containing 10 and 100 water molecules 
per interlaminar space, respectively. As can be seen, the 
time step shortens with increasing the system size, since 
energy fluctuations enlarge. This is why HMC is not very 
efficient for systems counting on a large number of mov- 
able sites. This is not our case, since few ions and water 
molecules are the only contributors to energy fluctua- 
tions. This makes HMC a reasonable choice. In fact, the 
time steps we are obtaining are similar to those usually 
employed for typical MD calculations^iiiiSSi^. 



For sampling in the NPzzT ensemble, after a trial 
change of particles' positions, a box change is attempted 
in such a way that the stress normal to the surface of 
the clays, Pzzi is kept constant. For this purpose, box 
fluctuations are allowed only in the z-direction and the 
probability for accepting the new box configuration is 
given by 

P = min{l,exp[-/3(AW + Pzz^V- 7V/3-iln(K/^o))]} 

(7) 

Here, /\U is the change in the potential energy, AV is 
the volume change, N is the total number of molecules, 
and Vn and Vo are the new and old box volumes, 
respectivelji^. 

For sampling in an open ensemble, the possibility of 
insertions and deletions of water molecules has to be con- 
sidered. Water insertions and deletions were performed 
by Rosenbluth samplingi2i2i. The ijlPzzT ensemble^ was 
used to obtain the equilibrium states when the system is 
in contact with a reservoir at certain temperature, pres- 
sure and water chemical potential. For this purpose, the 
algorithm must sample the probability density of finding 
the system in a particular configuration, i. e. 

.r y^exp{-/3[Z^(s^)-MiV+P^T/]} 
A/mP..t oc A3N^\ 

Hence, particle movements, insertions, deletions, and box 
changes must be done as in a typical NVT, ^VT and 
NPzzT sampling24. In particular, after trying a change 
of particles' positions, we performed tries of inserting- 
deleting water molecules. This is done by randomly call- 
ing both possible trials in such a way that calls are equally 
probable. Since accepting insertions or deletions are rare, 
we repeat this step 10 times or until any insertion or 
deletion is accepted. In case of refusing the 10 insertion- 
deletion trials, we performed a box trial move'^^. In this 
way, the system rapidly evolves to an equilibrium state. 
For some conditions, however, two free energy local min- 
ima appear, which are accessed by handling initial con- 
ditions. 



III. RESULTS 

Results are presented in two subsections. These are 
Sampling in the NPzzT ensemble and Sampling in the 
fiPzzT ensemble. Each part presents the results for 
ground level conditions, i. e., T— 298 K and P= 1 atm; 
and for 4 km of burial depth, i. e., T= 394 K and P= 
600 atm, assuming average gradients of 30 K/km and 150 
atm/km. 

A. Sampling in the NPzzT ensemble 

Let us first focus on the swelling behavior of the K- 
montmorillonite hydrates under ground level conditions. 
This is shown as a plot of the interlaminar distance as a 
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FIG. 1: Interlaminar distance as a function of the number of 
water molecules per interlaminar space. 



function of the number of water molecules obtained by 
NPzzT simulations. This curve starts from a dehydrated 
state having a interlaminar distance of 10.4 A, which is 
somewhat larger than the experimental value of 10.1 A^. 
For a small number of water molecules, the system reach 
interlaminar distances above 11 A, which slightly increase 
with increasing the number of water molecules produc- 
ing a plateau. This plateau shows interlaminar distances 
in the range 11.3 - 12.3 A, which are clearly lower than 
those produced by Na-montmorillonita^. This might 
be seen as something unexpected, since potassium ions 
are much larger than sodium ions. As we comment 
ahead, this is due to the different interlayer structures 
that sodium and potassium ions generate. For 60 water 
molecules per interlaminar space, there is a jump from a 
single water layer to a double water layer, which yields 
an interlaminar distance of 13.9 A. For larger amounts 
of water the system increases almost linearly, producing 
a very small step when jumping from a double to a triple 
water layer structure. 

Our findings are similar to those reported by Boek 
et al^. We included them in figure to make easy the 
comparison. It is seen that trends are practically equal. 
This indicates that differences in models and methods 
are not very important. Nevertheless, they always obtain 
slightly smaller interlaminar distances for a given number 
of water molecules. This difference is always lower than 
2%, except for the dehydrated state and the points close 
to 100 water molecules, where differences close to 5% are 
observed. Since in our case the jump from a single to 
a double layer is obtained for a lower amount of inter- 
laminar water, the points that correspond to 55-60 water 
molecules also show a larger difference. Anyway, differ- 
ences seem reasonable taking into account the differences 
in pair potential definitions, box sizes, and methodolo- 
gies. 

The structure of the interlaminar space for the system 
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FIG. 2: Oxygen, hydrogen and potassium density profiles of 
the interlaminar space. The water amount was fixed to 40 
molecules per interlaminar space. 



containing 40 water molecules per interlaminar space is 
shown in figure [21 A high and narrow oxygen peak at 
the interlaminar midplane, three hydrogen peaks, one 
coinciding with the oxygen peak and the other two sym- 
metrically situated at both sides, and a single potas- 
sium ion peak, also at the interlaminar midplane, are ob- 
served. Most of these results agree with others previously 
reportedi^iiSiii. It should be noted that this structure 
contrasts with that one of the Na-montmorillonite single 
layer hydrate. In this case, sodium ions are distributed 
at the sides of the oxygen peaks, closer to the clay sheets, 
some of them being strongly attached to the clay surface 
and so, forming inner-sphere surface complexes^. This 
makes water molecules to cluster in two layers joined at 
the interlaminar midplane, producing an oxygen double 
peak^^i^. Naturally, this effect widens the interlaminar 
space, despite of the smaller size of sodium ion. 

The radial distribution functions and coordination 
numbers for K-0 sites were also studied. Main g(r) peaks 
are situated at 2.77, 2.90, and 2.83 A for water, clay, and 
total oxygen sites, respectively. Although they are close 
to the one observed for bulk potassium water solution, 
a contraction is seen for the K-0 water distance as a 
consequence of confinement. On the other hand, the rel- 
atively large K-0 clay distance is due to the distribution 
of the potassium ions around the midplane of the inter- 
layer space. The coordination numbers for the first shell 
of oxygen atoms are 5.1, 5.0, and 10.1 for water, clay, 
and total oxygen sites, respectively. The clay and total 
coordination numbers are somewhat inflated due to the 
smaller separation of the 0-0 sites of the clay. This ex- 
plains why the total coordination number is larger than 
the one found for bulk. It should be noted the large con- 
tribution of the clay to the total coordination number. 
This suggests that potassium ions are interacting with 





FIG. 3: Snapshot of an equilibrated system having 40 water 
molecules. H sites are white, O are grey and K are black. The 
wireframe represents the clay structure. The topmost image 
is a side view and the lower image is the corresponding top 



both clay sheets at the same time. In addition, the large 
K-0 clay distance may be indicating that potassium ions 
are directly contributing to attract both layers toward 
the midplane, and hence holding the sheets close to one 
another. 

All the already pointed features are illustrated in fig- 
ure 13 For instance, it is observed that potassium ions 
are practically centered on the interlayer midplane with 
an approximate average of 5 water molecules surrounding 
each ion. It is also possible to see how ions tend to sepa- 
rate from each other, and that one of them is always close 
to a tetrahedral aluminum (these sites are not highlighted 
in the figure). Moreover, it is shown that there are no 
water molecules interposed between the potassium ions 
and the clay sheets. Hence, potassium ions are behaving 
as inner-sphere complexes, but simultaneously with both 
clay layers. To see this even clearer, figure 0] was built 
by rotating and zooming in figure |31 Here, a potassium 
ion and its inner water shell are shown. As can be seen, 
only 5 water molecules surround the ion that coordinates 
with two oxygen atoms from each clay. For this particu- 
lar case, are also seen average K-0 distances of 2.81 and 
3.04 A for water and clay, respectively. 

As mentioned before, for 60 water molecules a dou- 
ble layer hydrate is formed. This double layer hydrate 
becomes fully developed for 80 water molecules, where 




FIG. 4: Zoom in of a potassium ion and its coordination 
shell taken from figure |21 Only water molecules having K- 
O distances smaller than 3.7 A are shown. Distances in the 
figure are given in A. 
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FIG. 5: Oxygen, hydrogen and potassium density profiles of 
the interlaminar space. The water amount was fixed to 80 
molecules per interlaminar space. 



an interlayer distance close to 15.0 A is observed. For 
this system, the oxygen, hydrogen, and potassium pro- 
files are shown in figure [S] Here, the two oxygen peaks 
are a signature of the double water layer structure. For 
each oxygen peak a potassium peak and two hydrogen 
peaks are found. The potassium peak is found almost 
coinciding with the oxygen peak, but slightly displaced 
at the side closer to the clay sheet. This suggests the for- 
mation of inner-sphere complexes. On the other hand, 
a small hydrogen peak is found at the side closer to the 
clay sheet, and a larger one closer to the interlaminar 



6 




FIG. 6: Snapshot of an equilibrated system having 80 water 
molecules per interlaminar space. H sites are white, O are 
grey and K are black. The wireframe represents the clay 
structure. 



midplane. Hence, the third hydrogen peak found for the 
one layer case is missing here. 

For this configuration, main g(r) peaks are located at 
2.83, 2.75, and 2.79 A for water, clay, and total oxygen 
sites, respectively. Hence, the K-0 distance for water is 
not contracted anymore, but equal to the bulk water K- 
O distance. On the contrary, the K-0 distance for the 
clay decreased 0.17 A. Moreover, this distance is found 
to be almost constant for systems having more than 60 
water molecules, strongly suggesting that this is the nat- 
ural average K-0 distance for a potassium ion attached 
to the siloxane surface. Hence, the distance found for 
the single layer hydrate would be elongated as a result 
of the pressure the water molecules produce on the clay 
surfaces. The corresponding coordination numbers are 
6.8, 2.8, and 9.6 for water, clay, and total oxygen sites, 
respectively. As can be seen, the coordination number 
for the clay is much lower than the value found for the 
single layer hydrate. This means that potassium ions are 
coordinated to only one clay layer for the double layer hy- 
drate. Finally, and as expected, K-0 water coordination 
numbers increase with the number of water molecules, 
whereas K-0 clay and K-0 total coordination numbers 
decrease. A snapshot of this two layer hydrate is shown 
in figure ini We should mention that the structure of the 
inner-sphere complexes we found in this case are very 
similar to those already reported by Sposito et al}^ (not 
shown here). 

For burial conditions the interlaminar distance in- 
creases for a given number of water molecules per inter- 
laminar space. This is shown in figure [7| where it shows 
the interlaminar distance as a function of the number 
of water molecules for burial and for ground level con- 
ditions. This increment is more pronounced for large 
amounts of water. This is an expected behavior since 
water molecules occupy larger effective volumes at burial 
conditional. Consequently, profiles turn less sharp, i. e., 
peaks broad and shorten, as it is shown in figure |S1 
This was seen experimentally by Skipper et ai, although 
for a Na-montmorillonite systerr^. On the other hand, 
trends for both the ground and burial data are very sim- 
ilar, i. e., jumps are found at equal numbers of water 
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FIG. 7: Interlaminar distance as a function of the number 
of water molecules per interlaminar space for burial condi- 
tions. For comparison, the dotted line corresponds to ground 
conditions (from figure 0. 
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FIG. 8; Oxygen, hydrogen, and potassium density profiles 
of the interlaminar space. The water amount was fixed to 40 
molecules per interlaminar space for burial conditions. 



molecules. This last finding differs from the one found 
for Na-montmorillonite, where the single layer to double 
layer jump occurs for 60 or 50 water molecules, depend- 
ing on the burial depth^. 



B. Sampling in the ^PzzT ensemble 

Sampling in this ensemble allows the system to reach 
equilibrium with a reservoir whose temperature, pres- 
sure, and, in our case, water chemical potential are fixed. 
We employed the expression = /3/j,o -|-ln(p/po), where 
Po is the vapor pressure at equilibrium with liquid wa- 
ter whose chemical potential is /xqj E^nd p is the vapor 
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FIG. 9: Step evolution of the interlaminar distance and num- 
ber of water molecules obtained by fiP^^T sampling. In the 
lower plot, dashed lines correspond to each interlaminar space 
whereas the solid line is the average. Initial conditions were 
16 A and 60 water molecules per interlaminar space. Estab- 
lished conditions were T=298 K, P=l atm, and p/po=0.4. 
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FIG. 10: Interlaminar distance and number of water 
molecules per interlaminar space as a function of the vapor 
pressure for ground level conditions. Symbols □ and o corre- 
sponds to initial conditions of 10 water molecules and 11.5 A 
of interlaminar distance, and 60 water molecules and 16.0 A 
of interlaminar distance, respectively. 



pressure. For the TIP4P water model with T—298 K 
and P=l atm, we employed (3^o=-l7A. This value was 
obtained by NPT simulations of bulk TIP4P water— 
and using the Rosenbluth sampling method explained 
elsewhereii^. It is in good agreement with the value 
reported by Chavez-Paez et alJ^, and it is 2.5% larger 
than the one reported by Tambach et al^^ . For T=394 
K and P=600 atm, /3/io=-13.4 was obtained^^. 

The evolution of the interlaminar distance and num- 
ber of water molecules for ground level conditions and 
for pIpo=QA is shown in figurc|5| Here, the initial condi- 
tions were 16 A for the interlaminar distance and 60 wa- 
ter molecules per interlaminar space. It is observed that 
the system reaches approximately 65 water molecules and 
14.3 A of interlaminar distance at the initial simulation 
steps. This happens so fast since the system is initially 
very far from equilibrium. Immediately after, the system 
starts slowly loosing water molecules and decreasing its 
interlaminar distance on its way toward equilibrium. It 
is observed that, once the system reaches 13.5 A of in- 
terlaminar space, it quickly falls down to 12.2 A, loosing 
many water molecules during the process. This is a sig- 
nature of the transition from a double layer hydrate to 
a single layer hydrate. The amount of water that sig- 
nals the transition goes from 55 to 45 molecules. This 
agrees with the NPzzT results, where a single to a dou- 
ble layer transition is observed in the range of 40 - 60 
water molecules and 12.2 - 14.3 A. After the transition, 
it is observed that the system takes several steps to finally 
reach equilibrium at approximately 40000 steps. In this 
case, sampling was performed in the step range of 40000 
- 70000. For this particular run 12.0 A of interlaminar 
space and 39.3 water molecules were obtained. 

Many runs were performed in order to build up fig- 



ure 1101 Here it is seen the interlaminar distance and 
number of water molecules at equilibrium with reservoirs 
having different water vapor pressures and for ground 
level conditions. For zero vapor pressure, no matter what 
the established initial conditions, the system is forced to 
eliminate all of its water content, and so, the dehydrated 
state is yielded. This has 10.40 A of interlaminar space, 
as was found in the preceding section. For increasing 
the vapor pressure, the system starts to uptake water 
molecules from the reservoir producing larger interlam- 
inar distances. Nevertheless, the first water layer satu- 
rates for vapor pressures over 0.2po and so, a plateau is 
generated. Again, this single layer hydrate is produced 
for p/po < 0.4 for all initial conditions. 

For higher vapor pressures, however, at least two equi- 
librium states are yielded. That is, depending on the 
established initial conditions, the system may produce 
a double or a single layer hydrate. Hence, the topmost 
lines of both plots of figure ^| correspond to an initial 
configuration of 60 water molecules and the others to an 
initial configuration of 10 water molecules. These two 
equilibrium states are stable up to water vapor satu- 
ration, producing an open hysteresis cycle. The single 
layer state yields interlaminar distances ranging in 11.82 
- 12.06 A and amounts of water in the range of 30.1 - 41.8 
molecules. The double layer state produces interlaminar 
distances in the range 14.91 - 15.16 A and numbers of wa- 
ter molecules ranging in 76.7 - 80.6. These data contrast 
with those obtained for Na-montmorillonite2i. In this 
case, both initial conditions predict only a single layer hy- 
drate ioT p/po < 0.2. In addition, this single layer hydrate 
yields values of the interlaminar distance in the range 
of 12.34 - 12.73 A, having 34.6-47.5 water moleculesSi. 
Therefore, it becomes evident the K-montmorillonite ten- 
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dency to produce single layer hydrates even for relatively 
high vapor pressures, and the small interlaminar distance 
and amounts of water it yields. On the other hand, the 
double layer plateaus for potassium and sodium do not 
behave very differently for high vapor pressures. 

We should also compare our results with the experi- 
mental data obtained by Berend et alm^, Calvet'^®, and 
Sato et a/,^°. For this purpose, these data are also in- 
cluded in figure EH It should be pointed out that these 
data are direct measurements of interlaminar distances 
against the relative vapor pressure. This contrasts with 
the experimental data presented in figure ^ which were 
obtained by combining different experimental techniques 
and introducing some assumptions. Here, Berend et aZ.'s 
data were obtained as the evolution of the interlaminar 
distance for a dehydration process, by starting from a 
double or a single layer hydrate. Similarly, Calvet present 
his data by starting from a double or a single layer hy- 
drate, but for a hydration process. Sato et a/.'s data are 
also obtained by hydration. 

The best agreement between our simulations and ex- 
periments is produced when comparing with the data of 
Berend et alS^ We highlight the double layer hydrate for- 
mation for the range of relative vapor pressures of 0.6 - 
1.0; the single layer hydrate which keeps stable even for 
saturated vapor pressures; the very good agreement for 
the basal space for the single layer hydrate; and the rel- 
ative good agreement for the double layer hydrate inter- 
laminar distance, close to 15.5 A. Calvet's data also are 
in general agreement with our simulations. He observed 
a very stable single layer hydrate, with interlaminar dis- 
tances ranging in 11.8-12.5 A. Nevertheless, he obtained 
a very unstable double layer hydrate. Only for relative 
vapor pressures of 0.9 - 1.0 this hydrate was observed. Its 
interlaminar distance is in the range of 15.2-16.0 A. Sato 
et al.^s data also lead to similar interlaminar distances for 
the single and double layer hydrate. Finally, the single 
layer hydrate interlaminar distances also well agree with 
data reported by Brindlcy and Brown^. For 0.32, 0.52, 
and 0.79po of vapor pressure they reported 11.9, 11.9 and 
12.1 A, respectively. 

On the other hand, the most important discrepancy 
with experiments is that they obtained interlaminar dis- 
tances close to 10.0 A for the dehydrated state, which 
seems to be stable up to a relative vapor pressure rang- 
ing from 0.1 to 0.3. We think that if we were capable of 
reproducing the correct dehydrated distance, i. e. 10.0 in- 
stead of 10.4 A, this state would probably become stable 
for small vapor pressures, as they found. 

Experimentalists agree that for relative vapor pres- 
sures ranging in 0.0 - 0.4, a mixture of the dehydrated 
state and the single layer hydrate coexist, i. e. there is 
interstratification. Similarly, they observed the coexis- 
tence of double and single layer hydrates, for high rel- 
ative vapor pressures. In fact, Calvet refers to his own 
data as "apparent distances" ; Berend et al. conclude that 
all montmorillonites (they studied Li, Na, K, Rb and Cs- 
montmorillonites) form interstratified hydrates; and Sato 



TABLE II: g(r) main K-O peak positions, p, and first shell 
coordination numbers, n, for systems under different water 
vapor pressures, p/po- Subindexes w and c refer to water and 
clay, respectively. Peak positions are given in A. 





sin 


gle water layer 


double water layer 


p/po 


Pw 


Pc 


Tlyj 


ric 


Pm 


Pc 




0.1 


2.77 


2.83 


4.83 


5.32 








0.2 


2.77 


2.85 


4.96 


5.11 








0.3 


2.77 


2.88 


5.19 


4.94 








0.4 


2.77 


2.90 


5.17 


4.96 








0.6 


2.76 


2.89 


5.15 


4.93 


2.80 


2.75 


6.62 3.16 


0.8 


2.76 


2.88 


5.24 


4.85 


2.82 


2.77 


6.77 2.82 


1.0 


2.76 


2.90 


5.19 


4.89 


2.81 


2.75 


6.71 2.95 



et al. observed several "nonintegral basal reflections", 
which are interpreted as "random or segregated-type in- 
terstratification of collapsed and expanded layers". An 
example is the 13.81 A of interlaminar distance he ob- 
tained for the relative vapor pressure of 0.9. This point 
does not match either a single or a double hydrate in- 
terlaminar distance, as it is clearly seen in figure Uni It 
is important to mention that the general belief is that 
interstratification occurs due to chemical heterogeneities 
of the clay layers. It was already proven by simulations 
that changes on the positions of the clay substitutions 
produce different interlaminar distances (although dif- 
ferences are not very pronounced)!. Hence, these het- 
erogeneities surely lead to quasihomogeneous states, and 
probably, to interstratified ones. On the other hand, a 
perfect system like ours shows the double and the sin- 
gle layer hydrates to be stable for identical conditions. 
Hence, we do not see any reason for this not to occur in 
a real system. In other words, wc think that this is an- 
other source of interstratification, which arises just as a 
consequence of the inherent thermodynamics of the per- 
fect system. In fact, similar conclusions are deduced by 
Tambach et alw^. 

Table ^1 shows the main peak positions and the first 
shell coordination numbers for the systems equilibrated 
at different relative vapor pressures, and for a single and 
a double water layer configurations. This was done to 
prove that the shifts of the water and clay peaks men- 
tioned in the preceding section are indeed significant over 
a wide range of vapor pressures. As can be seen for the 
single layer hydrate, the K-0 water clay position is prac- 
tically constant and equal to the one found from the 
NPzzT sampling. Also the double layer is close to 
the value found in the previous section, i. e., similar to 
the K-0 bulk water position. On the other hand, pc in- 
creases for increasing the relative water pressure, reach- 
ing a plateau for p/po > 0.3. The plateau value is close 
to 2.89 A, which is also similar to the value found in 
the preceding section. The decrease of this distance with 
p/pq is a consequence of the decrease of the interlaminar 
distance with it. For the double layer configuration it is 
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FIG. 11: Interlaminar distance and number of water 
molecules per interlaminar space as a function of the vapor 
pressure for burial conditions. For comparison, the dotted 
line corresponds to ground conditions (from figure ITHll . 



also confirmed a value close to 2.76 A. Hence, the shifts 
of these peaks are relevant, and not just casual values 
obtained for some particular conditions. From table ITU it 
is seen that the values for the coordination numbers also 
agree with those presented in the previous section. 

To evaluate the effect of the burial depth on the K- 
montmorillonite swelling curves, figure ^] was built for 
T =394 K and P =600 atm. Additionally, results for 
ground level conditions were included to make the com- 
parison easy. As can be seen, for initial conditions of 10 
water molecules and a small interlaminar space, the inter- 
laminar distance is slightly smaller than those obtained 
for ground level conditions. This combines with the fact 
that water molecules occupy a larger effective volume for 
burial conditions and so, the number of water molecules 
decreases. In other words, the single layer hydrate of 
K-montmorillonite dehydrates under burial conditions. 
Furthermore, this single layer hydrate is stable even for 
a saturated water vapor pressure. Again, this contrasts 
with the results we obtained for Na-montmorillonite^^, 
where the single layer hydrate was found to be unsta- 
ble for large water activities. Moreover, for p/po = 1-0, 
the K-montniorillonite single layer hydrate yields 12.07 
A of interlaminar space and 40.4 water molecules, which 
seems to be far from the transition range of 45 to 55 
water molecules previously found. 

The double layer hydrate behaves differently. That 
is, it keeps constant its water content and it slightly in- 
creases its interlaminar distance. As found for ground 
level conditions, it collapses forming the single layer hy- 
drate for p/po < 0.6. On the other hand, for a satu- 
rated vapor pressure the system monotonically increases 
its amount of water as the simulation evolves. In this 
last case we stopped the simulations when reaching 180 
water molecules, and thus, we assumed that the system 
entered the osmotic regime. 



In summary, these results indicate that once formed 
the K-montmorillonite single water hydrate will not be 
destabilized neither at ground level conditions nor at 
burial depths. In addition, if the reservoir vapor pres- 
sure fells down 0.4po the single water hydrate will form. 



IV. DISCUSSION AND CONCLUSIONS 

Single and double layer hydrates of K-montmorillonite 
were studied by means of NP^zT and fiPzzT simulations. 
For that purpose a hybrid Monte Carlo scheme was em- 
ployed. Most results from the NPzzT sampling just con- 
firm those previously reported elsewhere^^'^^'^^-^'*-^^, sug- 
gesting that differences in models and methods are not 
very important. 

We found different clay-water-ion complexes for the 
single water hydrate. They consist of an average of 5 
water molecules surrounding a potassium ion that ad- 
ditionally coordinates with two oxygen atoms of each 
closest clay layer. These potassium ions are placed at 
the interlayer midplane and so, distances between them 
and the oxygens of the clay are slightly larger than those 
observed for the double layer hydrate. It should be men- 
tioned that these midplane ions were reported by pre- 
vious works as forming outer-sphere complexesi^ii^. In 
fact, Chang et al^ found that these midplane ions have 
a relatively large mobility, since potassium ions neither 
strongly coordinate to water molecules nor strongly inter- 
act with the clay surfaces. This is, indeed, characteristic 
of an outer-sphere complex. This explains why they call 
them in this way. Nonetheless, the inner-sphere definition 
given elsewherei^ says literally "The surface complex is 
inner-sphere if the cation is bound directly to a cluster of 
surface oxygen ions, with no water molecules interposed" . 
Accordingly, our midplane potassium ions form inner- 
sphere complexes but simultaneously with both clay lay- 
ers. Nevertheless, we expect them to have a diffusivity 
similar to that reported by Chang et ali^. This is a clear 
difference between this kind of inner-sphere complex and 
the typical inner-sphere complex, where ions are strongly 
coordinated to one of the surfaces^^. This simultaneous 
coordination of potassium ions with oxygen atoms of the 
two adjacent clay layers plus their relatively large coordi- 
nation distances, suggests to us that these complexes are 
attracting the clay layers toward the interlayer midplane, 
aiding to keep them together. 

Although the previous finding seems to explain the sta- 
bility of the single layer hydrate found experimentally, 
we should also reproduce computationally this behavior. 
Hence, /iP^^T simulations were carried out. This en- 
semble lets interchange water with a given reservoir and 
volume fluctuations. So, for any given water vapor pres- 
sure of the reservoir, an average amount of water and an 
average interlaminar distance are found. However, two 
different equilibrium states may be produced, pointing 
to the formation of free energy local minima^^, which are 
accessed by handling initial conditions. This produces 
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hysteresis loops2£. We found that for initial conditions 
close to the dehydrated state, a single water layer is al- 
ways obtained for any non zero vapor pressure, signatur- 
ing its high stability. The amount of water of this single 
layer slightly increases as the vapor pressure increases. 
This was obtained for ground level and burial conditions. 
The interlaminar distance for the ground level state is al- 
ways close to 12.0 A. The number of water molecules was 
found to be close to 38 for ground level conditions and 
about 36 for burial conditions. This is due to the larger 
effective volume the water molecules occupy at higher 
temperatures. 

On the other hand, the stability of the double layer hy- 
drate differs from that one of the single layer. For both 
conditions studied, it was seen that the double layer col- 
lapses to form the single layer for vapor pressures under 
0.4po- In addition, this double layer was found to be un- 
stable for burial depth and for saturated vapor pressures, 
producing a hydrated state in the osmotic regime. 

We should point out that these results agree with those 
of Boek et alm^. In their figure 3 is shown the potential 
energy of the interlayer water as a function of the num- 
ber of water molecules. It can be seen that the potassium 
curve is very different than the sodium and lithium ones 
only for small amounts of water (it shows much higher 
energy values, even well above the water bulk reference). 
This supports our finding of an extremely stable single 
water hydrate. Nevertheless, they claim that potassium 
is a good swelling inhibitor due to its ability to migrate 
and bind to the clay surfaces for large amounts of water. 
Hence, the negatively charged surface becomes screened, 
making its inherent repulsion less effective. We agree 
that this mechanism explains the relatively high stabil- 
ity of double layer hydrates, and even explains the stabil- 
ity of K-montmorillonite hydrates in the osmotic regime. 
Nevertheless, it cannot explain the remarkable stability 
of the single layer hydrate, since at least a double layer is 



needed to obtain the binding between ions and surfaces. 
This fact makes us think that midplane potassium ions 
are playing an important role in the stability of single 
layer hydrates. 

Despite of thinking that the potassium simultaneously 
binding to adjacent layers is the key to the single layer 
hydrate stability, we do not think it is the only factor. 
Lowest interaction energy for the pair K-TIP4P water is 
close to -20 kcal, which is higher than that for Na-TIP4P 
of approximately -25 kcal. This means that potassium 
ion does not strongly interact with water and therefore, 
it easily loses water from its coordination shell. In other 
words, the water chemical potential of the interlayer is 
not very negative. That is, water may be more comfort- 
able in bulk than in the interlayer, surrounding the ions. 
As mentioned by Sposito et al^, potassium ions show a 
kind of "hydrophobic character" , in the sense that they 
tend to interact with water molecules not only through 
their positive charge but also through solvent cage for- 
mation. Consequently, some of the water oxygen atoms 
may be easily exchanged by oxygen atoms from the clay 
surface. 

Finally, for drilling purposes and based on our results, 
it seems to not be enough to add potassium to the mud in 
order to guarantee the single layer water formation and, 
in this way, avoid swelling. One should simultaneously 
decrease the mud water activity to safely produce the 
single layer state. Once obtained, the clay will not swell 
at all. 
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